#include "checkVsTok.h"

void CheckVsT::read_Tok(std::string filename, int J2_t, int T2_t) {
  J2 = J2_t;
  T2 = T2_t;
  vector<int> data_i(9);
  double data_d;
  char u1[256];   // unuse
  std::string u2; // unuse
  std::ifstream file;
  file.open(filename);
  int size;
  file >> size;
  file.getline(u1, 256, '\n');
  file.getline(u1, 256, '\n');
  jaco_t.resize(size);
  for (int i = 0; i < size; i++) {
    file >> data_i[0] >> data_i[1] >> data_i[2] >> data_i[3] >> data_i[4] >>
        data_i[5] >> data_i[6] >> data_i[7] >> data_d >> data_i[8];
    for (int j = 0; j < 9; j++) {
      //  cout<<data_i[j]<<'\t';
    }
    // cout<<endl;
    jaco_t[i].J_2 = J2;
    jaco_t[i].T_2 = T2;
    jaco_t[i].T12_2 = data_i[1] * 2;
    jaco_t[i].n12 = data_i[2];
    jaco_t[i].l12 = data_i[3];
    jaco_t[i].S12_2 = data_i[4] * 2;
    jaco_t[i].J12_2 = data_i[5] * 2;
    jaco_t[i].n3 = data_i[6];
    jaco_t[i].l3 = data_i[7];
    jaco_t[i].E = data_i[8];
    jaco_t[i].j3_2 = data_d * 2;
  }
  file.close();
  if (-1 == 0) {
    for (int i = 0; i < size; i++) {
      cout << i + 1 << "\t";
      cout << jaco_t[i].T12_2 / 2 << "\t";
      cout << jaco_t[i].n12 << "\t";
      cout << jaco_t[i].l12 << "\t";
      cout << jaco_t[i].S12_2 / 2 << "\t";
      cout << jaco_t[i].J12_2 / 2 << "\t";
      cout << jaco_t[i].n3 << "\t";
      cout << jaco_t[i].l3 << "\t";
      cout << jaco_t[i].j3_2 / 2.0 << "\t";
      cout << jaco_t[i].E << "\t";
      cout << endl;
    }
  }
}

void CheckVsT::match(Jacobi &Jacobi_z) {
  // cout<<"$$$   "<<J2<<"\t"<<T2<<endl;
  int channel_num = 0;
  for (int i = 0; i < Jacobi_z.state.size(); i++) {
    if (Jacobi_z.state[i].size() > 0) {
      if (Jacobi_z.state[i][0].J_2 == J2 && Jacobi_z.state[i][0].T_2 == T2) {
        channel_num = i;
        // break;
      }
    }
  }
  cout << "match channel : " << channel_num << endl;
  channel_jac = Jacobi_z.FindC(J2, T2);
  channel = channel_num;
  int size_z = Jacobi_z.state[channel_num].size();
  int size_t = jaco_t.size();
  t2z.resize(size_t);
  z2t.resize(size_t);

  for (int i = 0; i < size_t; i++) {
    int flag = -1;
    t2z[i] = -1;
    for (int j = 0; j < size_z; j++) {
      if (Jacobi_z.state[channel_num][j].n12 != jaco_t[i].n12)
        continue;
      if (Jacobi_z.state[channel_num][j].n3 != jaco_t[i].n3)
        continue;
      if (Jacobi_z.state[channel_num][j].l12 != jaco_t[i].l12)
        continue;
      if (Jacobi_z.state[channel_num][j].l3 != jaco_t[i].l3)
        continue;
      if (Jacobi_z.state[channel_num][j].S12_2 != jaco_t[i].S12_2)
        continue;
      if (Jacobi_z.state[channel_num][j].J12_2 != jaco_t[i].J12_2)
        continue;
      if (Jacobi_z.state[channel_num][j].j3_2 != jaco_t[i].j3_2)
        continue;
      if (Jacobi_z.state[channel_num][j].T12_2 != jaco_t[i].T12_2)
        continue;
      if (Jacobi_z.state[channel_num][j].J_2 != jaco_t[i].J_2)
        continue;
      if (Jacobi_z.state[channel_num][j].T_2 != jaco_t[i].T_2)
        continue;
      if (Jacobi_z.state[channel_num][j].E != jaco_t[i].E)
        continue;
      t2z[i] = j;
      flag = 1;
    }
    if (flag == -1) {
      cout << "Can't Find this state :" << endl;
      cout << i + 1 << "\t";
      cout << jaco_t[i].T12_2 / 2 << "\t";
      cout << jaco_t[i].n12 << "\t";
      cout << jaco_t[i].l12 << "\t";
      cout << jaco_t[i].S12_2 / 2 << "\t";
      cout << jaco_t[i].J12_2 / 2 << "\t";
      cout << jaco_t[i].n3 << "\t";
      cout << jaco_t[i].l3 << "\t";
      cout << jaco_t[i].j3_2 / 2.0 << "\t";
      cout << jaco_t[i].E << "\t";
      cout << endl;
      // exit(0);
    }
  }
  // for(int i=0;i<size_t;i++)cout<<i<<"\t"<<t2z[i]<<endl;
}

void CheckVsT::read_Tok_TP(std::string filename, SingleP &SP, TripleP &TP) {
  std::ifstream file;
  file.open(filename);
  int i, na, la, nb, lb, nc, lc, J12, T12, E;
  double ja, jb, jc, J, T, MT;
  int size;
  int tz_a, tz_b, tz_c;
  int a_z, b_z, c_z;
  std::string unuse;
  char u1[256]; // unuse
  file >> size;
  file.getline(u1, 256, '\n');
  file.getline(u1, 256, '\n');
  TP_tok.resize(size);
  for (int k = 0; k < size; k++) {
    file >> i >> na >> la >> ja >> nb >> lb >> jb >> nc >> lc >> jc >> J12 >>
        J >> T12 >> T >> E;
    //cout<<i<<"\t"<<na<<"\t"<<la<<"\t"<<ja<<"\t"<<nb<<"\t"<<lb<<"\t"<<jb \
    //	<<"\t"<<nc<<"\t"<<lc<<"\t"<<jc<<"\t"<<J12<<"\t"<<J<<"\t"<<T12<<"\t" \
    //	<<T<<"\t"<<MT<<"\t"<<E<<endl;
    // cout<<k<<endl;
    int J_2 = J * 2;
    int T_2 = T * 2;
    int T12_2 = T12 * 2;
    int ch = TP.FindC(J_2, T_2);

    // int MT_2 = MT * 2;
    int ja_2 = ja * 2;
    int jb_2 = jb * 2;
    int jc_2 = jc * 2;
    // c_z=SP.FindSP(nc,lc,jc_2,tz_c);

    // a_z=SP.FindSP(na,la,ja_2,tz_a);
    // b_z=SP.FindSP(nb,lb,jb_2,tz_b);
    int index = -1;
    int size_tp = TP.state[ch].size();
    for (int j = 0; j < size; j++) {
      if (T12_2 == TP.state[ch][j].Tab_2 && J12 * 2 == TP.state[ch][j].Jab_2) {
        // if (E == TP.state[ch][j].E && MT_2 == TP.state[ch][j].T_2z) {
        if (E == TP.state[ch][j].E) {
          a_z = TP.state[ch][j].a;
          b_z = TP.state[ch][j].b;
          c_z = TP.state[ch][j].c;
          if (na == SP.state[a_z].n && nb == SP.state[b_z].n &&
              nc == SP.state[c_z].n) {
            if (la == SP.state[a_z].l && lb == SP.state[b_z].l &&
                lc == SP.state[c_z].l) {
              if (ja_2 == SP.state[a_z].j_2 && jb_2 == SP.state[b_z].j_2 &&
                  jc_2 == SP.state[c_z].j_2) {
                index = j;
                break;
              }
            }
          }
        }
      }
    }
    if (index == -1) {
      cout << "wrong at read_Tok_TP\t" << endl;
      cout << i << "\t" << na << "\t" << la << "\t" << ja << "\t" << nb << "\t"
           << lb << "\t" << jb << "\t" << nc << "\t" << lc << "\t" << jc << "\t"
           << J12 << "\t" << J << "\t" << T12 << "\t" << T << "\t" //<< MT << "\t"
           << E << endl;
      // cout<<a_z<<"\t"<<b_z<<"\t"<<c_z<<endl;
      cout << "Channel : " << ch << endl;
      exit(0);
    }
    TP_tok[k].channel = ch;
    TP_tok[k].index = index;
    TP_tok[k].T_2 = T_2;
    TP_tok[k].J_2 = J_2;
    /*
    cout<<k+1<<"\t";
    cout<<SP.state[a_z].n<<"\t"<<SP.state[a_z].l<<"\t"<<SP.state[a_z].j_2<<"\t";
    cout<<SP.state[b_z].n<<"\t"<<SP.state[b_z].l<<"\t"<<SP.state[b_z].j_2<<"\t";
    cout<<SP.state[c_z].n<<"\t"<<SP.state[c_z].l<<"\t"<<SP.state[c_z].j_2<<"\t";
    cout<<TP.state[ch][index].Jab_2/2.0<<"\t"<<TP.state[ch][index].J_2/2.0<<"\t";
    cout<<TP.state[ch][index].Tab_2/2<<"\t";
    cout<<TP.state[ch][index].T_2/2.0<<"\t"<<TP.state[ch][index].T_2z/2.0<<"\t";
    cout<<TP.state[ch][index].E<<endl;
    */
  }
  file.close();
}

void CheckVsT::read_Tok_T(std::string filename, Jacobi_NL &jacNL, int T_2) {

  /*
   check_t.read_Tok("data/Tok_Jac.dat",1,1);
   check_t.match(jaco_b);
   check_t.read_Tok_T("data/Tok_T.dat",jaco_nl);

   size=check_t.T_tok.size();
   for(int i=0;i<size;i++){
     int ch_TP=check_t.T_tok[i].ch_TP;
     int ch_NL=check_t.T_tok[i].ch_NL;
     int x=check_t.T_tok[i].index_TP;
     int y=check_t.T_tok[i].index_NL;
     int E=triple_t.state[ch_TP][x].E;
     int x_index = x - triple_t.channel[ch_TP].indexE_beg[E];
     int y_index = y - jaco_nl.channel[ch_NL].indexE_beg[E];
     //cout<<check_t.T_tok[i].i<<"\t"<<check_t.T_tok[i].p<<endl;
     // cout<<y<<"\t index jac "<<check_t.T_tok[i].index_jac<<"\t"<<x<<endl;
     //cout<<"good?\t"<<ch_TP<<"\t"<<E<<"\t"<<y_index<<"\t"<<x_index<<endl;
     double T=trans_t.Tval_vec[ch_TP][E][y_index][x_index];
     //cout<<check_t.T_tok[i].i<<"\t"<<check_t.T_tok[i].p<<"\t"<<T<<endl;
   }

   cout<<"test T
   "<<trans_t.cal4(jaco_nl.state_NL[1][7],triple_t.state[1][6])<<endl;
   cout<<"test T "<<trans_t.Tval_vec[1][1][7][2]<<endl;
   */

  std::ifstream file;
  file.open(filename);
  int i, p, N0, L0;
  double T;
  char u1[256]; // unuse
  // file.getline(u1, 256, '\n');
  T_tok.clear();
  // cout<<"good?"<<endl;
  while (file.peek() != EOF) {
    file >> i >> p >> N0 >> L0 >> T;
    if (file.eof())
      break;
    // cout<<i<<"\t"<<p<<endl;
    i = i - 1;
    p = p - 1;
    int index_jac = t2z[p];
    int ch_TP = TP_tok[i].channel;
    // cout<<"good?"<<endl;
    int ch_NL = jacNL.FindC(TP_tok[i].J_2, TP_tok[i].T_2);
    if (TP_tok[i].T_2 != T_2) {
      cout << "wrong at read TOk T for T_2" << endl;
      cout << i + 1 << "\t" << p << "\t" << N0 << "\t" << L0 << endl;
      cout << TP_tok[i].J_2 << "\t" << TP_tok[i].T_2 << "\t" << ch_NL << endl;
      cout << index_jac << "\t" << N0 << "\t" << L0 << endl;
      exit(0);
    }
    // cout<<"good?"<<endl;
    int size = jacNL.state_NL[ch_NL].size();
    int index_NL = -1;
    for (int k = 0; k < size; k++) {
      if (index_jac == jacNL.state_NL[ch_NL][k].index) {
        if (channel_jac == jacNL.state_NL[ch_NL][k].JT_channel) {
          if (N0 == jacNL.state_NL[ch_NL][k].N_0 &&
              L0 == jacNL.state_NL[ch_NL][k].L_0) {
            index_NL = k;
            break;
          }
        }
      }
    }
    if (index_NL == -1) {
      // cout<<"wrong at read TOk T"<<endl;
      // cout<<i+1<<"\t"<<p<<"\t"<<N0<<"\t"<<L0<<endl;
      // cout<<TP_tok[i].J_2<<"\t"<<TP_tok[i].T_2<<"\t"<<ch_NL<<endl;
      // cout<<index_jac<<"\t"<<N0<<"\t"<<L0<<endl;
      // exit(0);
    }
    T_state_Tok T_t;
    T_t.index_TP = TP_tok[i].index;
    T_t.ch_TP = ch_TP;
    T_t.ch_NL = ch_NL;
    T_t.index_NL = index_NL;
    T_t.index_jac = index_jac;
    T_t.N0 = N0;
    T_t.L0 = L0;
    T_t.T_tok = T;
    T_t.i = i + 1;
    T_t.p = p + 1;
    if (index_NL == -1) {
      T_t.flag = -1;
    } else {
      T_t.flag = 1;
    }
    //
    T_tok.push_back(T_t);
  }
  cout << "T_tok size = " << T_tok.size() << endl;
}

void CheckVsT::check_cont(Cont_term &cont_x, int J2_c, int T2_c) {

  int channel_c = channel_jac;
  cout << "Channel : " << channel_c << endl;
  cout << "J2 = " << J2_c << " \t T2 = " << T2_c << endl;

  for (int i = 0; i < this->t2z.size(); i++) {
    for (int j = 0; j < this->t2z.size(); j++) {
      int x = t2z[i];
      int y = t2z[j];
      double res = 0.0;
      if (x >= 0 && y >= 0) {
        res = cont_x.vec[channel_c][x][y];
      }
      if (std::abs(res) > 0.00001) {
        cout << i + 1 << "\t" << j + 1 << "\t" << x << "\t" << y << "\t" << res
             << endl;
      }
    }
  }
}

void CheckVsT::check_cont(Cont_term &cont_x, int J2_c, int T2_c,
                          string filename) {

  std::ofstream file;
  file.open(filename);

  int channel_c = channel_jac;
  file << "Channel : " << channel_c << endl;
  file << "J2 = " << J2_c << " \t T2 = " << T2_c << endl;

  for (int i = 0; i < this->t2z.size(); i++) {
    for (int j = 0; j < this->t2z.size(); j++) {
      int x = t2z[i];
      int y = t2z[j];
      double res = 0.0;
      if (x >= 0 && y >= 0) {
        res = cont_x.vec[channel_c][x][y];
      }
      if (std::abs(res) > 0.00001) {
        file << i + 1 << "\t" << j + 1 << "\t" << x << "\t" << y << "\t" << res
             << endl;
      }
    }
  }
  file.close();
}

void CheckVsT::check_ope(Ope_term &ope_x, int J2_c, int T2_c) {

  int channel_c = channel_jac;
  cout << "Channel : " << channel_c << endl;
  cout << "J2 = " << J2_c << " \t T2 = " << T2_c << endl;

  for (int i = 0; i < this->t2z.size(); i++) {
    for (int j = 0; j < this->t2z.size(); j++) {
      int x = t2z[i];
      int y = t2z[j];
      double res = 0.0;
      if (x >= 0 && y >= 0) {
        res = ope_x.vec[channel_c][x][y];
      }
      if (std::abs(res) > 0.00001) {
        cout << i + 1 << "\t" << j + 1 << "\t" << x << "\t" << y << "\t" << res
             << endl;
      }
    }
  }
}

void CheckVsT::check_ope(Ope_term &ope_x, int J2_c, int T2_c, string filename) {

  std::ofstream file;
  file.open(filename);

  int channel_c = channel_jac;
  file << "Channel : " << channel_c << endl;
  file << "J2 = " << J2_c << " \t T2 = " << T2_c << endl;

  for (int i = 0; i < this->t2z.size(); i++) {
    for (int j = 0; j < this->t2z.size(); j++) {
      int x = t2z[i];
      int y = t2z[j];
      double res = 0.0;
      if (x >= 0 && y >= 0) {
        res = ope_x.vec[channel_c][x][y];
      }
      if (std::abs(res) > 0.00001) {
        file << i + 1 << "\t" << j + 1 << "\t" << x << "\t" << y << "\t" << res
             << endl;
      }
    }
  }
  file.close();
}

void CheckVsT::check_JacA_vec(Package &pack, string filename) {
  std::ofstream file_A;
  file_A.open(filename);
  file_A << "!!! Cation the print is 3 times of vec result" << endl;
  // int size_ch=JacA_vec.size();
  // for(int ch=0; ch<size_ch; ch++){
  // int size=JacA_vec[ch].size();
  file_A << " Channel : " << channel_jac << endl;
  file_A << "J2 = " << J2 << " \t T2 = " << T2 << endl;
  for (int i = 0; i < t2z.size(); i++) {
    for (int j = 0; j < t2z.size(); j++) {
      int x = t2z[i];
      int y = t2z[j];
      double jacA_res = 0.0;
      if (x >= 0 && y >= 0) {
        jacA_res = pack.JacA_vec[channel_jac][x][y] * 3;
      }
      if (std::abs(jacA_res) > 0.00001)
        // file<<i<<"\t"<<j<<"\t"<<JacA_vec[ch][i][j]*3<<endl;
        file_A << i + 1 << "\t" << j + 1 << "\t" << x << "\t" << y << "\t"
               << jacA_res << endl;
    }
  }

  file_A.close();
}

void CheckVsT::check_Dcoeff(Jacobi &jaco_b, Antisym &antisym_t,
                            string filename) {
  std::ofstream file_D;
  file_D.open(filename);
  file_D << " Channel : " << channel_jac << endl;
  file_D << "J2 = " << J2 << " \t T2 = " << T2 << endl;
  file_D << "i\tj\ta\tb\t E_i\t E_j\t D_coefficients" << endl;
  for (int i = 0; i < t2z.size(); i++) {
    for (int j = 0; j < t2z.size(); j++) {
      int x = t2z[i];
      int y = t2z[j];
      int E_x = jaco_b.state[channel_jac][x].E;
      int E_y = jaco_b.state[channel_jac][y].E;
      int x_index = -1;
      int y_index = -1;
      double ant_res = 0.0;
      if (x >= 0 && y >= 0) {
        if (E_x == E_y) {
          x_index = x - jaco_b.channel[channel_jac].indexE_beg[E_x];
          y_index = y - jaco_b.channel[channel_jac].indexE_beg[E_y];
          ant_res = antisym_t.vec_D[channel_jac][E_x][x_index][y_index];
        }
      }
      if (std::abs(ant_res) > 0.000001)
        file_D << i + 1 << "\t" << j + 1 << "\t" << x << "\t" << y << "\t"
               << E_x << "\t" << E_y << "\t" << ant_res << "\t" << x_index
               << "\t" << y_index << endl;
    }
  }
  for (int i = 0; i < antisym_t.eig_val[channel_jac][5].size(); i++) {
    // cout<<i<<"\t"<<antisym_t.eig_val[channel_c][5][i]<<endl;
  }

  file_D.close();
}

void CheckVsT::check_Tcoeff(TripleP &triple_t, Jacobi_NL &jaco_nl,
                            TransCoeff &trans_t, string filename) {
  std::ofstream file_T;
  file_T.open(filename);
  file_T << "J2 = " << J2 << "\t T2 = " << T2 << endl;
  int size = T_tok.size();
  for (int i = 0; i < size; i++) {
    // cout<<"i : "<<i<<endl;
    // cout<<"\t"<<check_t.T_tok[i].i<<"\t"<<check_t.T_tok[i].p<<"\t"<<endl;
    // if(check_t.T_tok[i].flag == -1)continue;
    int ch_TP = T_tok[i].ch_TP;
    int ch_NL = T_tok[i].ch_NL;
    int x = T_tok[i].index_TP;
    int y = T_tok[i].index_NL;
    //   cout<<"good 1?"<<endl;

    // cout<<ch_TP<<"\t"<<ch_NL<<endl;
    // cout<<x<<"\t\t"<<y<<endl;

    int E = triple_t.state[ch_TP][x].E;
    // if (E != jaco_nl.state_NL[ch_NL][y].E)continue;

    double T = 0.0;
    if (x >= 0 && y >= 0) {
      int x_index = x - triple_t.channel[ch_TP].indexE_beg[E];
      int y_index = y - jaco_nl.channel[ch_NL].indexE_beg[E];
      // cout<<"good 2?"<<endl;

      // cout<<check_t.T_tok[i].i<<"\t"<<check_t.T_tok[i].p<<"\t"<<check_t.T_tok[i].N0;
      // cout<<"\t"<<check_t.T_tok[i].L0<<endl;

      // cout<<ch_TP<<"\t"<<ch_NL<<"\t"<<E<<"\t"<<x<<"\t"<<y<<endl;
      // cout<<"y_index="<<y_index<<"\t x_index="<<x_index<<" \t E = "<<E<<endl;

      T = trans_t.Tval_vec[ch_TP][E][y_index][x_index];
    }
    // cout<<"good 3?"<<endl;
    double diff = T - T_tok[i].T_tok;
    file_T << T_tok[i].i << "\t" << T_tok[i].p << "\t";
    file_T << T_tok[i].N0 << "\t" << T_tok[i].L0 << "\t";
    file_T << ch_TP << "\t" << x << "\t" << ch_NL << "\t" << y << "\t";
    file_T << T_tok[i].T_tok << "\t" << T << "\t";
    if (std::abs(diff) > 0.00001) {
      cout << T_tok[i].i << "\t" << T_tok[i].p << "\t";
      cout << T_tok[i].N0 << "\t" << T_tok[i].L0 << "\t";
      cout << ch_TP << "\t" << x << "\t" << ch_NL << "\t" << y << "\t";
      cout << T_tok[i].T_tok << "\t" << T << "\t";
      cout << diff << endl;
      file_T << "\t" << diff;
    }
    file_T << endl;
    // cout<<y<<"\t index jac "<<check_t.T_tok[i].index_jac<<"\t"<<x<<endl;
    // cout<<"good?\t"<<ch_TP<<"\t"<<E<<"\t"<<y_index<<"\t"<<x_index<<endl;
    // cout<<check_t.T_tok[i].i<<"\t"<<check_t.T_tok[i].p<<"\t"<<T<<endl;
  }
  file_T.close();
}

void CheckVsT::check_pack(Package &pack, string filename) {
  std::ofstream file;
  file.open(filename);
  int size = TP_tok.size();
  file<<"i\t j\t x\t y\t channel\t ME(MeV)"<<endl;
  for (int k = 0; k < size; k++) {
    int ch = TP_tok[k].channel;
    for (int kk = 0; kk < size; kk++) {
      int x = TP_tok[k].index;
      int y = TP_tok[kk].index;
      // cout<<x<<"\t"<<y<<endl;
      int ch_k = TP_tok[k].channel;
      int ch_kk = TP_tok[kk].channel;
      if (ch_k != ch_kk)
        continue;
      if (std::abs(pack.vec[ch][x][y]) > 0.000001) {
        file << k + 1 << "\t" << kk + 1 << "\t" << x << "\t" << y << "\t" << ch
             << "\t" << pack.vec[ch][x][y] << endl;
      }
    }
  }
  file.close();
}

void CheckVsT::read_Tok_tpe(std::string filename, int J2_t, int T2_t) {
  std::ifstream file;
  file.open(filename);
  std::string unuse;
  char u1[256]; // unuse
  file.getline(u1, 256, '\n');
  file.getline(u1, 256, '\n');
  file.getline(u1, 256, '\n');
  file.getline(u1, 256, '\n');
  Tpe_struct tpe_t;
  while (!file.eof()) {
    int p, q;
    double C1, C3, C4;
    file >> p >> q >> C1 >> C3 >> C4;
    tpe_t.p = p;
    tpe_t.q = q;
    tpe_t.C1 = C1;
    tpe_t.C3 = C3;
    tpe_t.C4 = C4;
    cout << p << "\t" << q << endl;
    Tpe_tok_vec.push_back(tpe_t);
  }
  file.close();
}

void CheckVsT::check_tpe(Tpe_term &tpe_x, int J2_c, int T2_c) {

  int channel_c = channel_jac;
  cout << "Channel : " << channel_c << endl;
  cout << "J2 = " << J2_c << " \t T2 = " << T2_c << endl;

  int size_Tpe = Tpe_tok_vec.size();
  cout << "p\t q\t C1_tok \t C1_zhuo\t C1_diff\t C3_tok \t C3_zhuo\t "
          "C3_diff\t "
          "C4_tok \t C4_zhuo\t C4_diff\t"
       << endl;
  for (int t = 0; t < size_Tpe; t++) {
    int i = Tpe_tok_vec[t].p - 1;
    int j = Tpe_tok_vec[t].q - 1;

    int y = t2z[i];
    int x = t2z[j];
    double C1, C3, C4;
    if (x >= 0 && y >= 0) {
      C1 = tpe_x.vec1[channel_c][x][y];
      C3 = tpe_x.vec3[channel_c][x][y];
      C4 = tpe_x.vec4[channel_c][x][y];
    }
    double C1_t = Tpe_tok_vec[t].C1;
    double C3_t = Tpe_tok_vec[t].C3;
    double C4_t = Tpe_tok_vec[t].C4;
    double diff_c1 = std::abs(C1 - C1_t);
    double diff_c3 = std::abs(C3 - C3_t);
    double diff_c4 = std::abs(C4 - C4_t);
    if (std::abs(diff_c1) < 0.000001)
      diff_c1 = 0.0;
    if (std::abs(diff_c3) < 0.000001)
      diff_c3 = 0.0;
    if (std::abs(diff_c4) < 0.000001)
      diff_c4 = 0.0;
    // if (std::abs(C1 + C3 + C4) > 0.00001) {
    cout << i + 1 << "\t" << j + 1 << "\t" //<< y << "\t" << x << "\t"
         << std::scientific << std::setprecision(2) << C1_t << "\t" << C1
         << "\t" << diff_c1 << "\t" << C3_t << "\t" << C3 << "\t" << diff_c3
         << "\t" << C4_t << "\t" << C4 << "\t" << diff_c4 << endl;
    //}
  }
}
